home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / SeqPups / apps / clustalw.src / random.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  1.5 KB  |  79 lines  |  [TEXT/R*ch]
  1. /*
  2. *    
  3. *    Rand.c
  4. *    
  5. *    -    linear and additive congruential random number generators
  6. *        (see R. Sedgewick, Algorithms, Chapter 35)
  7. *
  8. *    Implementation: R. Fuchs, EMBL Data Library, 1991
  9. *    
  10. */
  11. #include <stdio.h>
  12.  
  13. unsigned long linrand(unsigned long r);
  14. unsigned long addrand(unsigned long r);
  15. void addrandinit(unsigned long s);
  16.  
  17. static unsigned long mult(unsigned long p,unsigned long q);
  18.  
  19.  
  20. #define m1    10000
  21. #define m    100000000
  22.  
  23. /* linear congruential method
  24. *    
  25. *    linrand() returns an unsigned long random number in the range 0 to r-1
  26. */
  27.  
  28.  
  29. unsigned long linrand(unsigned long r)
  30. {
  31.     static unsigned long a=1234567;
  32.     
  33.     a = (mult(a,31415821)+1) % m;
  34.     return( ( (a / m1) * r) / m1 );
  35. }
  36.  
  37. static unsigned long mult(unsigned long p, unsigned long q)
  38. {
  39.     unsigned long p1,p0,q1,q0;
  40.     
  41.     p1 = p/m1; p0 = p % m1;
  42.     q1 = q/m1; q0 = q % m1;
  43.     return((((p0*q1 + p1*q0) % m1) * m1 + p0*q0) % m);
  44. }
  45.  
  46.  
  47. /* additive congruential method
  48. *    
  49. *    addrand() returns an unsigned long random number in the range 0 to r-1
  50. *    The random number generator is initialized by addrandinit()
  51. */
  52.  
  53. static unsigned long j;
  54. static unsigned long a[55];
  55.  
  56. unsigned long addrand(unsigned long r)
  57. {
  58. int x,y;
  59. /*        fprintf(stdout,"\n j = %d",j);  */
  60.     j = (j + 1) % 55;
  61. /*        fprintf(stdout,"\n j = %d",j);  */
  62.     x = (j+23)%55;
  63.     y = (j+54)%55;
  64.     a[j] = (a[x] + a[y]) % m;
  65. /*    a[j] = (a[(j+23)%55] + a[(j+54)%55]) % m;  */
  66. /*        fprintf(stdout,"\n a[j] = %d",a[j]);     */
  67.     return( ((a[j] / m1) * r) / m1 );
  68. }
  69.  
  70. void addrandinit(unsigned long s)
  71. {
  72.     a[0] = s;
  73.     j = 0;
  74.     do {
  75.         ++j;
  76.         a[j] = (mult(31,a[j-1]) + 1) % m;
  77.     } while (j<54);
  78. }
  79.